function [stdz_th,stdst_th,irfz,irfst,acz0_th,acz1_th,vdecz_th,vdecst_th,vdechz,vdechst,vcvz_th]=...
    irfandthmom(Gone,RRone,SDXone,ZZ,nperirf,npervdec,addzvec,flag_unit)
% =========================================================================
% MODE_MOMENTS 
% Compute theoretical and simulation moments for a model solution passed
% If flag_unit==1 then compute the IRF to a UNIT impulse 
% else, compute to a 1 std impulse 
% Alejandro Justiniano (c)
% =========================================================================
ny=size(Gone,1);nx=size(SDXone,1);nz=size(ZZ,1); 
if nargin <  9
    flag_unit=0; 
else
    if flag_unit~=1 && flag_unit~=0
        flag_unit=0; 
    end
end
flag_addz=0;
if nargin >= 8 && ~isempty(addzvec)
    flag_addz=1;
    if length(addzvec)~=nz
        error('ADDZVEC must be NZ x 1')
    end
    indic_addz=find(addzvec==1); 
end
% =============================
%% 1.B IRF Baseline Model 
% =============================
irfz=zeros(nz,nperirf+1,nx); 
irfst=zeros(ny,nperirf+1,nx); 
for ii=1:nx 
    if flag_unit==0
        resp=SDXone(ii,ii)*RRone(:,ii);
    else
        resp=RRone(:,ii);
        disp('Unit Impulse'); 
    end
    irfst(:,1,ii)=resp; 
    irfz(:,1,ii)=ZZ*resp; 
    for jj=2:(nperirf+1); 
        resp=Gone*resp;
        irfst(:,jj,ii)=resp;
        irfz(:,jj,ii)=ZZ*resp;
    end 
    if flag_addz==1
        % Corrected July 14 2008 
        irfz(indic_addz,:,ii)=cumsum(squeeze(irfz(indic_addz,:,ii)),2);
    end
end 
% ===========================================================
% Variance Deconposition at different horizons as well as 
% Theoretical Moments for Check 
% ===========================================================
[vdecz_th,vcvz_th,stdz_th,acz0_th,acz1_th,junk,junk,vdechz,vdechst,vdecst_th,stdst_th]=vardecom_mode(Gone,RRone,SDXone,ZZ,npervdec); 
%stdz_th=sqrt(diag(vcvz_th)); 